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ABSTRACT 

We compute the two-, three- and four-point correlation functions from the Wilkinson Microwave 
Anisotropy Probe ( WMAP) first-year data, and compare these to a Monte Carlo ensemble of 5000 
realizations, based on the best-fit WMAP running-index spectrum of Gaussian fluctuations. The 
analysis is carried out in three steps, covering small (< 72'), intermediate (< 5°) and large scales (up 
to 180°). On the largest scales our results are consistent with the previously reported hemisphere 
power asymmetries: the northern ecliptic hemisphere is practically devoid of large scale fluctuations, 
while the southern hemisphere show relatively strong fluctuations. We also detect excess correlations 
in W-band difference maps as compared to the detailed noise simulations produced by the WMAP 
team, possibly indicative of unknown systematics. While unlikely to affect any temperature based 
results, this effect could potentially be important for the upcoming polarization data. On intermediate 
angular scales we find hints of a similar anisotropic distribution of power as seen on the very largest 
scales, but not to the same extent. In general the model is accepted on these scales. Finally, the same 
is also true on the smallest scales probed in this paper. 

Subject headings: cosmic microwave background — cosmology: observations — methods: statistical 
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1. INTRODUCTION 

In recent months, a large number of analyses focus- 
ing on non-Gaussianity in the Wilkinson Microwave 
Anisotropy Probe {WMAP; Bennett et al. 2003a) data 
have clai med signific a nt de tec tions of non-Gaussian 
featiir es (iConi et alJ 120041: Ide Oliveira- Costa et all 
20041 lEriksen et alJ I20p4albl: IHansen et all l2004alht 
Larson fc Wandelt' '20041: iMcEwen et alJ 120041: iPard 
2004: Vielva et al...2004) 7 If any one of these detections 
can be shown to be of cosmological origin, currently 
accepted models based on Gaussianity and isotropy 
would have to be revised. Gaining proper understanding 
of their nature is therefore essential for further progress. 

Sources of non-Gaussian (or anisotropic) signal may 
be categorized into three general classes. First, most 
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non-cosmological foregrounds are highly non-Gaussian, 
and are all likely to introduce a non-Gaussian signal 
into the maps to some extent. In fact, unless some 
particular detection is explicitly demonstrated to be 
frequency-independent, it must usually be assumed to 
be foreground-induced. Second, systematics may intro- 
duce non-Gaussian signals into the data. An example 
of this is correlated noise, that results in stripes along 
the scanning path of the experiment. Finally, the most 
intriguing possibility is that the primordial density field 
itself could be non-Gaussian, e.g., through the existence 
of topological defects or non-equilibrium infiation. 

In the current paper, we subject the WMAP data to 
an analysis based on real-space A'^-point correlation func- 
tions. While harmonic-space methods often are preferred 
over real-space methods for studying primordial fluctu- 
ations, real-space methods may have an advantage with 
respect to systematics and foregrounds, since such ef- 
fects are usually localized in real space. It is therefore 
important to analyze the data in both spaces in order 
to highlight different features. For instance, by consid- 
ering difference maps between independent differencing 
assemblies (abbreviated DA; see Hinshaw et al. [2003] for 
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details on the terminology), that ideally should contain 
no CMB signal, we detect excess correlations in the data 
that are not accounted for in detailed simulations of the 
WMAP pipeline, and by partitioning the sky into small 
regions, we find hints of residual foregrounds near the 
galactic plane. 

The al gorithms used in th is paper were devel- 
oped by lEriksen et alJ ll2004cl^. and apphe d to the 
first-year WMAP data bv lEriksen et all l|2004a|) . 
Other A^-point correlation function analyses of the 
first-year WMA P data in clude those presente d by 
CTaztafiaga et al.l|2 003). Ga ztaiiaga fc Wagd l|2n0.'^ . and 
Land fc Magueiiol (|2004) . 



2.1. 



2. DEFINITIONS 

The statistics of interest in this paper are the TV-point 
correlation functions (here restricted to two-, three- and 
four-point functions), and we measure these functions 
both for the observed data and for an ensemble of simu- 
lated realizations with controlled properties. A statis- 
tic is then employed to quantitatively measure the agree- 
ment between the data and the model. 

An A^-point correlation function is by definition the 
average product of N temperatures, measured in a fixed 
relative orientation on the sky, 
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where the unit vectors ni, . . . , fiAr span an A'^-point poly- 
gon on the sky. By assuming statistical isotropy, the Ap- 
point functions are only functions of the shape and size of 
the A^-point polygon, and not on its particular position 
or orientation on the sky. Hence, the smallest number of 
parameters that uniquely determines the shape and size 
of the A^-point polygon is 2A^ — 3. 

The Af-point correlation functions are estimated by 
simple product averages. 



^2N 



-3) 



"N 



rpi rp 
■ It^ - ■ ■ 1^ 



N 



(2) 



"N 



where the sums are taken over all sets of A^ pixels ful- 
filling the geometric requirements set by 9i, . . . , 92N-3- 
The pixel weights, Wi, may be independently chosen for 
each pixel in order to reduce, e.g., noise or border ef- 
fects. Here they represent masking, by being set to 1 for 
included pixels and to for excluded pixels. 

The main difficulty with computing A"-point functions 
is their computational scaling. The number of inde- 
pendent pixel combinations scales as 0{N^^^), and for 
each combination of A^ pixels, 2A^ — 3 angular distances 
must be computed to uniquely determine the properties 
of the corresponding polygon. Computing the full Ap- 
point function for N > 2 and A'pi^ > lO'^ is therefore 
computationally challenging. 

However, it is not necessary to include all possible 
A'^-point configurations in order to produce interesting 
results. E.g., one may focus only on small angular 
scales, or on configurations with some special symme- 
try propertie s . By using the methods described by 
lEriksen et al.l l)2004cD . the computational expense then 
becomes tractable, since no CPU time is spent on ex- 
cluded configurations. In this paper several such subsets 
are computed, covering three distinct ranges of scales, 
namely small (up to 1?2), intermediate (up to 5°) and 
large scales (the full range between 0° and 180°). 



The statistic 

In this paper, a simple test is chosen to quantify 
the degree of agreement between the simulations and the 
observations, where as usual is defined by 
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iCNii) - {CN{i)))M^\CN{j) - {Cn{j))). 



(3) 

Here C^ii) is the A^-point correlation function for config- 
uration'* number i, (CAr(i)) is the corresponding average 
from the Monte Carlo ensemble, and 



^^^"-^ k=l 

(4) 

is the covariance matrix. 

This statistic is optimized for studying Gaussian dis- 
tributed data. Unfortunately, the A^-point correlation 
functions (and in particular even-ordered ones) are gen- 
erally strongly non-Gaussian (and asymmetrically) dis- 
tributed, and this leads to an uneven weighting of the 
two tails by the statistic. In order to remedy this 
weakness, the empirical distributi on of each configura- 
tion is transformed by the relation ijEriksen et alJl2004a(l 



Rank of observed map 
Total number of maps + 1 



1 



r^* dt. 



(5) 



The numerator of the left hand side is the number of 
realizations with lower value than the current map, and 
the denominator is the total number of realizations plus 
1. The addition of 1 is necessary to obtain symmetric 
values of s around 0, and to avoid that the realization 
with the lowest value is assigned an infinite confidence 
level. Note that if the data were in fact Gaussian dis- 
tributed, equation 13 is an identity operation in the limit 
of an infinite number of simulations. The statistic is 
then computed from the transformed data, rather than 
from the original correlation functions. 

The quoted significance level is given in terms of the 
fraction of simulations with a lower x^ value than the 
observed map. Thus, a value more extreme than either 
0.025 or 0.975 indicates that the model is rejected at the 
2 (7 level. 

In order to eliminate any procedural difference between 
the simulations and the observed maps, we include the 
observed map itself in the estimation of the covariance 
matrix. While this should have no impact on the result if 
the covariance matrix is properly converged, it is a very 
useful safe-guard against such issues. 

A singular value decomposition (SVD) is used to com- 
pute the inverse covariance matrix, and conservatively all 
modes with a condition number smaller than 10^^ are set 
to zero. However, this limit is only reached in the small- 
scale analysis, in which different neighboring configura- 
tions are very strongly correlated, and the covariance ma- 
trix converges more slowly than for the intermediate- and 
large-scale functions. 

Finally, the four-point correlation function is treated 
differently than the two- and three-point functions, in 

The terms 'configuration' and 'bin' are used interchangeably 
in this paper. 



3 




Fig. 1. — The low-resolution co-added WMAP map is made 
by smoothing each of the eight cosmologically interesting bands 
to a common FWHM = 140' Gaussian beam, and subsequently 
co-adding these using inverse-noise weights. Finally, best-fit 
monopole, dipole, and quadrupole moments were removed from 
the high-latitude region. 



that its power spectrum dependence is reduced by utiliz- 
ing the foUowing relationship: if a random field is Gaus- 
sian, then the ensemble average of the four-point function 
may be written in terms of the two-point function (see, 
e.g., Adler 1981), 

(TiTaTgTi) ={TiT2){T3Ti) + {T^T^i) {T2Ti) (6) 

+ {TiTi){T2T3). (7) 

We therefore subtract the quantity on the right-hand side 
from the observed four-point function, to obtain a re- 
duced four-point function. In what follows, all results 
for the four-point function refer to this reduced function. 

3. PREPARATION OF THE DATA 

The first-year WMAP data may be downloaded from 
LAMBDA^ . Most of the analyses described in the follow- 
ing sections are carried out for both the raw maps a nd 
the template-corrected versions ijBennett et al.ll2003bj) . 

We define our model for the simulations as the sum 
of a CMB component and a noise component. The sig- 
nal component is based on the best-fit WMAP power 
spectrum with a running index, including multipolc 
components with I = 2, . . . , 1024, filtered through the 
HEALPix*^ (Gorski, Hivon, & Wandelt 1999) pixel win- 
dow functions and channel-dependent beam windows. 
While there is some controversy about the evidence for 
a running index, we have found that this spectrum pro- 
vides a better fit to the data at the very low-^ range 
of the spectrum, and therefore a better fit in terms of 
iV-point correlation functions that are sensitive to large- 
scale structures. However, this difference is only notice- 
able for the two-point function; the three-point and re- 
duced four-point functions are only mildly dependent on 
the assumed power spectrum, thus our results should be 
independent thereof. Finally, the aim's are assumed to 
be Gaussian. 

The noise is assumed to be uncorrelated and Gaussian, 
with rms levels given for each pixel of each channel by 
the WMAP team (Bennett et al. 2003a). This noise is 
added pixel by pixel and channel by channel to the CMB 
signal realizations. 

^ http://lambda.gsfc.nasa.gov 

® http:/ /www. eso.org/science/healpix 



We study both individual frequency maps and a 
CO- added version that includes all eight bands. The 
frequency maps are generated by straight averaging 
over bands using equal weights, whereas the co-added 
map is weight ed with inverse noise variance weights 
(jHinshaw et al.L2003a) . 

The analysis is carried out in two steps: First we study 
the large-scale fluctuations on the full sky^ by degrading 
the maps from Ngidc — 512 to Ngidc = 64, Ngidc being the 
HEALPix resolution parameter (Gorski et al. 1999). We 
then study the small and intermediate scales by parti- 
tioning the full-resolution sky into disks of 10° radius (in 
two different configurations), and compute the correla- 
tion functions on each disk separately. Full-sky functions 
for these scales are estimated by averaging over all disks. 

The degradation process may be written on the follow- 
ing algorithmic form: 

1. Compute the spherical harmonic components, a^m 
from the full resolution A^sido = 512 map. 

2. Deconvolve with the original WMAP beam and 
pixel windows (i.e., multiplication in harmonic 
space). 

3. Convolve with a 140' FWHM Gaussian beam and 
Aside — 64 pixel windows. 

4. Compute the A"sidc = 64 map using the filtered 

This process is carried out for each channel separately 
before any co-addition is done. The downgraded, co- 
added WMAP map is shown in Figure ^ 

Since all structures in the high-resolution maps are 
smoothed out in the degrading process, the foreground 
exclusion mask must also be extended correspondingly. 
This is done by setting all excluded pixels in the original 
mask to and all included pixels to 1, then convolving 
this map with a Gaussian beam of the desired FWHM, 
and finally excluding all pixels with a value smaller than 
0.99 in the smoothed mask. 

This degraded mask will only be used on smoothed, 
low- resolution maps, and the KpO mask is therefore 
used with point sources not excluded as our base mask. 
This could in principle introduce a non-Gaussian sig- 
nal into our maps, but in practice point sources con- 
tribute negligi ble pow er at scales larger than a few de- 
grees ( Hinsha w et al.ll200 3a) . 

Finally, for the low-resolution analysis, we remove 
the monopole, dipole, and quadrupole modes from each 
map separately, with parameters computed from the 
high-latitude regions of the sky only (defined by the 
extended KpO mask). The reason for removing the 
quadrupole is that thi s particular mode m ay have an 
anomalously low value l)Bennett et al.ll2'003a|) , but is cer- 
tainly contaminated by residual foregrounds afte r tem- 
plate s ubtraction (Slosar & Seljak 2004; Hanse n et alJ 
20043 lEriksen etT a l. 2004d). Since real-space correla- 
tion function are inherently more sensitive to the low-^ 

Whenever we refer to a 'full-sky' analysis, we mean that data 
from both hemispheres are included, except for those pixels excised 
to avoid contamination from the galactic plane and point-sources, 
where appropriate. 
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Fig. 2. — The large-scale, low-resolution correlation functions computed from from the co-added WMAP map. The solid dots show the 
results from the observed data, the solid line and the gray bands show the median and the 1, 2, and 3 a confidence regions, respectively, 
computed from 5000 simulations. The full-sky, northern, and southern galactic hemisphere results are plotted in left, middle, and right 
columns, while rows show the two-point, the equilateral three-point, and the rhombic four-point functions. Note in particular the extremely 
featureless correlation functions compute d on the northern hem isphere, indicating little large-scale structure in this region. Similar plots 
for the ecliptic hemispheres are shown by lEriksen et al.l I2004al) . 



modes, this well-known effect could mask other interest- 
ing features. 

In the case of the small- and intermediate-scale analy- 
ses, we estimate the correlation functions on independent 
disks of 10° radius. In order to reduce the correlation be- 
tween neighboring disks, we therefore choose to remove 
all multipoles^ with < ^ < 18, by generalizing the usual 
method of removing the low-£ components to higher mul- 
tipoles. 

4. LARGE-SCALE ANALYSIS 

The first analysis focuses on the very largest scales, by 
computing the A'^-point correlation functions from de- 
graded maps, as described above. The functions are uni- 

* The particular ^max = 18 was chosen to correspond roughly to 
the disk radius of 10° . 



formly binned with 1° bin size, and the two-point func- 
tion is computed over the full range betwe en and 180°. 
For the higher-order functions we follow lEriksen et alJ 
(,2002) and compute the pseudo-coUapsed and the equi- 
lateral three-point functions, and the 1+3-point and the 
rhombic four-point functions. 

The definition of "pseudo-coUapsed" is slightly mod- 
ified compared to the one described by lEriksen et alJ 
1,2002) . In this paper "pseudo-coUapsed" indicates that 
the length of the collapsed edge falls within the second 
bin, and not that only neighbori ng pixels are included 
ijGaztanaga et aH2003t lEriksen et al. 200A^. This mod- 
ification eliminates the need for treating each configura- 
tion as a special case, and is thus purely implementation- 
ally motivated. 

The results from these measurements are shown in Fig- 
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Fig. 3. — Distributions of the values computed from the ensemble for the full sky three-point function; the plots show the results for 
the northern ecliptic hemisphere (left), the southern ecliptic hemisphere (middle), and the full sky (right). The value corresponding to the 
foreground corrected co-added map is marked with a solid line, while the raw co-added map is denoted with a dashed line. 



ureHlfor the co-added map, for a few selected functions. 
A complete summary of the large-scale measurements are 
given in Tabled for both individual channels and for the 
co-added map, and for two different masks. 

Considering first the full-sky two-point function, we see 
that this function demonstrates an almost complete lack 
of structure above 60°, and its overall shape is very flat 
as pointed out by several authors (e.g., Bennett et al. 
2003a). However, here it is important to remember that 
the quadrupole was removed prior to the computation 
of the correlation functions, and therefore the two-point 
function does not appear quite as anomalous as that seen 
in many other plots^. 

Next, the full-sky three-point function shows similar 
tendencies, as it lies inside the 1 a confidence region al- 
most over the full range of scales. Finally, the four-point 
function is quite low at small angles, and very close to 
zero at large angles. Thus, all three full-sky correlation 
functions point toward the same conclusion - there is 
little large-scale power in the WMAP data. 

Several analyses have presented evidence for a sig- 
nificant asymmetry between the norther n and south- 
ern ecliptic (and galactic) hemispheres jEriksen et al.l 
|2p04ab; Hansen et al. 2004a b; Park 200J), and there- 
fore we choose to estimate the various functions from 
these regions separately. Similar patterns are indeed 
found also in these cases: the northern hemisphere cor- 
relation functions all show a striking lack of fluctuations, 
whereas the southern hemisphere functions show good 
agreement with the confidence bands computed from the 
Gaussian simulations. This difference translates into a 
clear difference of numbers, as seen in Tabled The 
northern hemisphere results for the higher-order func- 
tions all lie in the bottom few percent range, while the 
corresponding numbers for the southern hemisphere are 
generally higher than 80%. 

^ Although the two-point correlation function is the Legendre 
transform of the power spectrum, it does not necessarily follow that 
the observed two-point function agrees with an ensemble average 
based on a power spectrum fitted to the data: the best-fit power 
spectrum is largely determined by the small scale information (high 
^'s) in the data, whereas the two-point function is very sensitive 
to the largest scales (low Ps). The two functions thus provide 
complementary pictures of the data, highlighting different features. 



The statistic may actually serve as a general mea- 
sure of the overall fluctuation level of the higher-order 
functions, since they both have vanishing mean (we use 
the reduced, not the complete, four-point function in 
these analyses; for the same reason, this does not work 
for the two-point function). One possible statistic for the 
degree of power asymmetry between two complimentary 
hemisphere is therefore simply the ratio of the two in- 
dividual x^'s- This quantity is computed for both the 
simulations and the observed data, and the fraction of 
simulations with a smaller ratio is listed in the third row 
of each section in Table ^ 

We see that this ratio is extreme at the few percent 
level for the three-point function, and at less than one 
percent for the four-point function, for the ecliptic hemi- 
spheres. Further, it is not particularly sensitive to fre- 
quency or galactic cut. In fact, the numbers are slightly 
stronger for the conservative |6| > 30° mask than for the 
KpO mask in four-point function case. Both these results 
argue strongly against a foreground based explanation. 

In order to quantify the two-point function asymmetry, 
we adopt a slightly modified ve rsion of the S'-statistic, as 
defined bv lSpergel et alJ l|2003D 

j [C2{e)fdcose. (8) 

Note that we choose to include the full range of available 
angles, while iSpergel et alJ l)2003D chose to exclude an- 
gles smaller than 60°. Excluding the smaller angles does 
increase the nominal significance of this statistic when 
applied to the WMAP data, but it also makes the inter- 
pretation of the final results less clear, since the cut-off 
scale is arbitrarily chosen. 

In general, the S'-statistic has similar properties to a 
X^ statistic that only includes diagonal terms, but it has 
a distinct advantage over the latter in the case of the 
two-point function: while both power deficits and ex- 
cesses lead to a large x^ (rendering this statistic useless 
for probing asymmetry), the opposite is true for the S'- 
statistic. Power deficit yields a low S- value, while power 
excess yields a high S-value. In other words, this statistic 
may serve the same purpose for the two-point function 
as the x^ statistic does for the higher-order functions. 

The results from this analysis are shown in the second 
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TABLE 1 

Large scale Af-poiNT correlation function and 5-statistic results 

Q-band V-band VF-band Co-added 

Region KpO \b\ > 30° KpO |b| > 30° KpO |b| > 30° KpO |6| > 30° 



Two-point function; statistic 



Full sky 


0.725 


0.558 


0.491 


0.606 


0.519 


0.429 


0.574 


0.451 


Northern galactic 
Southern galactic 


0.495 
0.903 


0.605 
0.861 


0.721 
0.792 


0.732 
0.865 


0.682 
0.740 


0.574 
0.629 


0.772 
0.865 


0.578 
0.803 


Northern ecliptic 
Southern ecliptic 


0.439 
0.272 


0.508 
0.729 


0.608 
0.215 


0.582 
0.676 


0.218 
0.125 


0.536 
0.572 


0.538 
0.216 


0.469 
0.598 


Two-point function; 5-statistic results 


Pull sky 


0.094 


0.102 


0.084 


0.121 


0.068 


0.131 


0.085 


0.107 


Northern galactic 
Southern galactic 
Ratio of S-values 


0.026 
0.739 
0.033 


0.439 
0.433 
0.239 


0.019 
0.776 
0.025 


0.379 
0.555 
0.168 


0.036 
0.726 
0.034 


0.417 
0.567 
0.203 


0.022 
0.749 
0.029 


0.419 
0.488 
0.214 


Northern ecliptic 
Southern ecliptic 
Ratio of 5-values 


0.016 
0.745 
0.045 


0.015 
0.693 
0.075 


0.017 
0.791 
0.038 


0.012 
0.784 
0.052 


0.012 
0.743 
0.048 


0.013 
0.794 
0.050 


0.015 
0.759 
0.044 


0.014 
0.735 
0.063 






Three-point function; x^ 


statistic 








Pull sky 


0.314 


0.711 


0.267 


0.636 


0.154 


0.616 


0.284 


0.641 


Northern galactic 
Southern galactic 
Ratio of x^'s 


0.041 
0.825 
0.030 


0.051 
0.796 
0.046 


0.036 
0.819 
0.027 


0.050 
0.847 
0.033 


0.034 
0.819 
0.027 


0.061 
0.774 
0.057 


0.031 
0.822 
0.026 


0.060 
0.821 
0.044 


Northern ecliptic 
Southern ecliptic 
Ratio of x'^'s 


0.047 
0.831 
0.031 


0.046 
0.792 
0.040 


0.023 
0.861 
0.014 


0.033 
0.820 
0.031 


0.014 
0.871 
0.012 


0.038 
0.829 
0.031 


0.034 
0.840 
0.023 


0.041 
0.793 
0.039 






Pour-point function; x^ 


statistic 








Pull sky 


0.484 


0.518 


0.491 


0.480 


0.474 


0.472 


0.468 


0.508 


Northern galactic 
Southern galactic 
Ratio of x^'s 


0.089 
0.884 
0.030 


0.073 
0.920 
0.022 


0.069 
0.905 
0.020 


0.053 
0.930 
0.014 


0.086 
0.878 
0.031 


0.051 
0.914 
0.017 


0.077 
0.888 
0.025 


0.061 
0.923 
0.018 


Northern ecliptic 
Southern ecliptic 
Ratio of x^'s 


0.070 
0.852 
0.030 


0.020 
0.927 
0.004 


0.054 
0.873 
0.021 


0.010 
0.942 
0.001 


0.050 
0.846 
0.025 


0.012 
0.931 
0.002 


0.058 
0.857 
0.025 


0.014 
0.932 
0.002 



Note. — Results from tests of the large-scale correlation functions. The numbers 
indicate the fraction of simulations with a value lower than for the respective WMAP 
map. 



section of Tabled Overall, they are consistent with the 
X^-based three- and four-point function results, with the 
single exception of the galactic \b\ > 30° measurements, 
which do not show any signs of asymmetry. However, 
in this case the two-point function is quite poorly con- 
strained at the largest angles because of the limited sky 
coverage, and sample variance dominates the statistic. 

In Figure |21 the histograms of the values for the 
co-added simulated ensemble are plotted together with 
the observed WMAP values (both for the foreground- 
corrected and the raw maps). In this Figure it is well 
worth noting the effect of foregrounds, namely that the 
increases if foregrounds are present. This is both 
an intuitive and an important result: it is intuitive be- 
cause the statistic basically measures the amount of 
deviations from the average function, and for a func- 
tion with vanishing mean such as the three-point func- 
tion, it therefore quantifies the overall level of fluctua- 
tions. By adding a statistically independent component 



to the maps (residual foregrounds in our setting), more 
fluctuations are introduced into the three-point func- 
tion. This observation is therefore also important, since 
it implies that residual foregrounds are unlikely to ex- 
plain the northern hemisphere anomaly - sub-optimal 
foreground templates would introduce large-scale fluctu- 
ations, rather than suppress them. This also suggests 
that one could use the statistic as defined above to 
fit for the template amplitudes, a possibility that will be 
explored further in a future publication. 

All in all, the results presented in this section seem 
to disfavor a foreground-based explanation for the large- 
scale power asymmetry. The variation from band to band 
is very small indeed, and similar signs of asymmetry can 
be seen in any one of the frequencies. Further, there is 
no clear dependence on the particular sky cut. 

4.1. Analysis of difference maps 
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Fig. 4. — Results from the difference map analysis. The solid dots show the results computed from the observed maps, while the solid 
line and the gray bands show the median and the la confidence region, respectively, computed from the 110 simulations produced by the 
WMAP team, including all known systematic effects. The dashed lines indicate the 1 cr confidence region assuming uncorrelated noise, 
modulated by N^haip) only. 



Next, we study the noise properties of the WMAP 
data. SpecificaUy, the two-point correlation functions 
are computed from all possible difference maps within 
each frequency channel, and we compare these to the 
functions computed from corresponding correlated noise 
maps^° produced by the WMAP team. All data have 
been processed in the same way as for the large-scale 
analysis; the maps are downgraded to a common reso- 

Available at |http: / /lambda. gsfc.nasa.gov| 



lution of 140', then the difference maps are formed, and 
finally, the best-fit monopole, dipole and quadrupole mo- 
ments are removed from the extended KpO region. 

In Figure 0] the results from this analysis are shown, 
comparing the WMAP data (solid dots) to the simu- 
lations, that includes all known systematics. The gray 
bands indicate the 1 a confidence bands computed from 
the 110 available simulations, and the dashed lines show 
the 1 a regions assuming uncorrelated (but inhomoge- 
neous) white noise, computed from 1000 simulations. 
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In particular two features stand out in these plots: 
First, there is a strong peak at 6h — 141°, the effective 
horn separation angle of the WMAP satellite. Second, 
there is a clear rise toward high values at small angles, 
which is a real-space manifestation of correlated noise. 
Of course, neither of these effects are unexpected, since 
they are also present in the si mulations, and they ar e 
both discussed at some length bv lHinshaw et al.1 1)20031:)|) . 

However, there are a few surprises to be found in the 
M^-band plots. Specifically, a very strong signal may be 
seen in the Wl-WA map. To the extent that the confi- 
dence regions can be approximated by Gaussians, we see 
that the peak at 9h extends to more than 4 a compared 
to the simulations, and the overall fluctuation levels are 
clearly stronger than what is seen in the simulations. A 
similar pattern is also seen in the W1~W2 map, but with 
a slightly smaller amplitude. 

Comparing the confidence regions estimated from cor- 
related and white noise simulations, we see that the main 
difference is more large-scale curvature in the correlated 
noise bands. This is consistent with the power spec- 
trum view, where correlated noise is found to have the 
strongest impact at low £'s. This again translates into a 
two-point function with a shape resembling that of the 
signal-dominated functions shown in Figure |2| This ef- 
fect is particularly evident in the maps which involve 
the W4 differencing assembly, which is known to have a 
significantly higher knee frequency than the other DA's 
(Ijarosik et al. 2003). 

We now quantify the agreement between the observa- 
tions and the model by means of a statistic, but we 
do not attempt to include the correlation structure in 
this case because of the limited number of simulations. 
Rather, we define a simplified statistic on the following 
form, 

i=l ^ ' 

Here A^bin is the number of bins in the correlation func- 
tion, and (C2{i)) and cr^(i) are the average and variance, 
respectively, of bin number i computed from the simula- 
tions. Unsurprisingly, this statistic strongly rejects the 
model for the WI-W2 and Wl-WA combinations, as 
none of the 110 simulations have a higher x^j^g value 
than the observation, or even close to it. For the re- 
maining six combinations, the ratios of simulations with 
a lower Xdiag comfortably in the range between 

0.28 and 0.94. 

It is difficult to make firm conclusions about the origin 
of these structures based on this simple analysis alone, 
but it is evident that the noise simulations do not fully 
capture the nature of the data. On the other hand, it 
is also very unlikely that this effect has any significant 
impact on the cosmological results from the first-year 
WMAP data release, given its relatively small amplitude. 
It may be important with respect to the second-year po- 
larization data. 

A si milar detection was reported bv lFosalba fc Szapudil 
l)2004|) . They found that the noise contribution in the 
WMAP data may have been underestimated by 8-15% in 
the original analysis. However, it is difficult to establish a 
direct connection between these results, considering that 
their results are most significant at high €'s, while our 



analysis is explicitly restricted to low £'s. 

5. SMALL- AND INTERMEDIATE-SCALE ANALYSIS 

In order to probe smaller scales, subsets of the A^-point 
functions are now from the full-resolution co-added map. 
This analysis is facilitated by partitioning the sky into 
non-overlapping disks of 10° radius, each including be- 
tween 15 000 and 25 000 pixels (the number varies be- 
cause of the Kp2 mask). Two different sets of disks (de- 
noted A and B) are used in the following analysis, their 
union covers a total of 81% of the sky, and about 60% 
each. The two sets contain respectively 87 and 81 disks. 

The reasons for dividing the sky into patches are two- 
fold: First, the computational cost soon becomes difficult 
to handle for data sets with more than about 150000 pix- 
els. Since the algorithms scale as a relatively high power 
of A^pix, it is much cheaper to divide the full region into 
patches. Second, and equally important, we want to be 
able to localize interesting effects in pixel space. In par- 
ticular, we seek to study the effects discussed in ^2] fur- 
ther, and one convenient way of doing this is by analyz- 
in g the sky in patches. A similar analysis was carried out 
bv lHansen et al.l l)2004bD , using a power spectrum based 
statistic. 

The disk sets are created as follows: In each set, the 
disks are laid out on rings of constant latitude, with as 
many disks on each ring as there is space for without 
overlap (the polar rings of set B are exceptions to this 
rule), and random initial longitude. Then the WMAP 
Kp2 mask is applied and we keep only those "disks" (at 
this point some have a rather peculiar geometry) with 
more than 15 000 accepted pixels. The reason for prefer- 
ring the more liberal Kp2 mask over KpO is that we also 
want to study the effect of foregrounds in this analysis. 

The defining difference between set A and B is given 
by the latitudes on which the disks are centered. In set 
A, the disks are laid out on latitudes given hy 6 — k- 20°, 
k = 0, . . . , 9, while the disks in set B are centered on 
61 = (fc + 1/2) • 20°, fc = 0, 8. This difference implies 
that set A has two rings of disks that touches the galactic 
plane, while the center ring of disk set B is completely 
discarded. It is therefore reasonable to assume that disk 
set A is more affected by foregrounds than disk set B, 
as will be confirmed later. The two disk sets are shown 
in the two top panels of Figure |31 superimposed on the 
co-added WMAP map. 

5.1. Intermediate- scale analysis 

First, we consider the A^-point correlation functions on 
intermediate scales, here defined as scales smaller than 
5-10°. Each function is binned with T.2 bin size, and the 
two-point function is computed up to 10°, for a total of 
83 bins. The three-point function is computed over all 
isosceles triangles for which the baseline is the longest 
edge, but no longer than 5°. Note that this set includes 
the equilateral triangle and three points on a line as spe- 
cial cases. Finally, the four-point function is computed 
over the same set of configurations, but with a fourth 
point added by reflecting the third point about the base- 
line. The total number of independent configurations is 
about 460. Note that since there are many more isosce- 
les triangles with 5° baseline than with 1°, the vast ma- 
jority of these configurations span scales from 3° to 5°. 
Consequently, the following analysis is dominated by 
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Fig. 5 . Results from the intermediate scale correlation function analysis. The top panels show the layout of each disk set, and the other three rows show the x 

results. The colors indicate the confidence level at which the disk is accepted, computed according to equation^^ Thus, dark blue indicates a very low x value, green a 
value around the median, and dark red a very high value. 
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intermediate scales rather than smaU scales, even though 
a few small-scale configurations are included. 

The results from this intermediate scale analysis are 
plotted in Figure where the colors indicate the confi- 
dence level at which each disk is accepted, as computed 
by equation However, extreme limits of —2.5 a and 
2.5(7 are enforced, because we only have a limited num- 
ber of simulations available. Here it is worth recalling 
that we removed all power with £ < 18 from the maps, 
and neighboring disks are therefore nearly uncorrelated. 
Distributions of the confidence level distributions are 
shown in Figure El 

We see that the two-point function is to a very good ap- 
proximation accepted by the test. There are no visibly 
connected patches of similar values, and the distribution 
of confidence values appears to be typical compared to 
the simulations. The three-point functions are more sus- 
picious, especially when considering the pattern seen in 
disk set B. In this case, two large, connected patches 
of low values are visible on the northern hemisphere, 
while the south-east quadrant appears to have quite large 
values. In other words, the asymmetry pattern found 
in the large-scale functions is apparent even in this plot. 
For disk set A, these features are less clear, but still con- 
sistent with disk set B; the northern hemisphere on aver- 
age have quite low values (or little fiuctuations), while 
the south-east quadrant has quite high values. This is 
particularly evident if one disregards all disks touching 
(i.e., are partially cut by) the galactic plane, they are 
more likely to be affected by residue foreground contam- 
ination. 

The four-point functions are less decisive, and the gen- 
eral agreement with the Gaussian model seems to be 
quite good. No particular features are seen in these cases. 

We now estimate the full sky correlation functions by 
averaging over all the individual disk correlation func- 
tions, weighting each sub- function by the number of pixel 
combinations, Nc, 

Here Cn is the full-sky iV-point correlation function, Cj^ 
is the j'th disk correlation function, and i represents the 
geometric configuration under consideration. 

This function is computed from each disk set individu- 
ally and over the union of the two sets. While the latter 
function obviously has the advantage of larger sky cov- 
erage, it also weights configurations that are completely 
contained in the intersection of two disks twice. On the 
other hand, the number of such common configurations 
is fairly small, at least in this intermediate-scale analysis. 

In Figure [71 the full-sky, intermediate-scale two-point, 
equilateral three-point and rhombic four-point functions 
are plotted, as computed from equation 1101 The corre- 
sponding x^ results are shown in Table [21 We see that 
the agreement between observations and simulations is 
in all cases very good, both in terms of overall shape 
and amount of small-scale fluctuations. Note, however, 
that these figures only show a small subset of the config- 
urations included in the full analysis; while there are 41 
three-point configurations with a base line of 5°, there are 
only two with a 10' base line. Thus, the right hand sides 
of the three- and four-point function plots are weighted 



much more strongly than the left hand side in the x^ 
analysis. Generally speaking, visualizing higher-order 
functions is difficult because of their high dimensionality. 

The results from the corresponding x^ analysis are 
shown in Table [3 Here we see that the results for the 
foreground corrected map all lie comfortably between 
0.05 and 0.95, and no sign of discrepancy is found. For 
the raw maps, the x^ numbers are generally somewhat 
high, but not disconcertingly so. The agreement with 
the assumed model on intermediate scales appears to be 
satisfactory on intermediate scales, as far as A^-point cor- 
relation functions are concerned. 

A similar analysis, including disks in the galactic or 
ecliptic hemispheres only, was also performed, but it did 
not find any clear discrepancy in either case. Thus, 
the asymmetric pattern seen in Figure [SI in the three- 
point function in set B does not correspond directly to a 
dipole type distribution. Of course, we could subdivide 
the sky further according to the observed patterns, but 
this would strongly dilute the final probabilities, since 
we then define our test a-posteriori. All in all, the in- 
termediate scale correlation functions accepts the model, 
although hints of hemisphere asymmetry may be seen by 
eye in the three-point function. 

5.2. Small-scale analysis 

The analysis from the previous section is now repeated, 
but this time including all three-point configurations 
with a longest edge shorter than or equal to 72' (about 
220 different configurations), and twice as many four- 
point configurations. The four-point configurations are 
defined in terms of the three-point configurations, by let- 
ting the fourth point either be mirrored or rotated about 
the base line, that once again is defined to be the longest 
edge of the triangle. 

The results from this disk-based analysis are shown 
in Figure [HI and distributions of the corresponding x^ 
values are plotted in Figure [HI By eye, the three-point 
function results in Figure [SI appear to be in quite good 
agreement with the model, and no clear anomalies stand 
out. This impression is confirmed both by the full-sky 
X^ numbers, as well as by the plots showing the disk x^ 
distributions, except for the fact that there are quite a 
large number of disks in the 1-1.5 a range in disk set B. 

However, the four-point function plot for disk set A 
shows a more interesting effect; at least 7 out of the 21 
disks touching the galactic sky cut in disk set A have 
a fairly high x^ value, and there are no disks with low 
X^ values. This is most likely an indication of residual 
foregrounds near the galactic plane, a conclusion which 
becomes eve n more plausib l e when studying figure 11 in 
the paper bv lBennett et al.l l|2008bt) . In these plots, clear 
residuals are seen outside the Kp2 mask, particularly in 
the Q-band map. 

We may quantify the significance of this effect by com- 
puting a new disk-averaged correlation function. This 
time we include only those 21 disks in the two near- 
galactic rows (A34~54), and the corresponding results 
are shown in Table 01 for both intermediate and small 
scales. 

The four-point function results have a combined sig- 
nificance at 2 CT for the intermediate scales, and almost 
3 cr for the small scales. In fact, for the small scales even 
the three-point function has a x^ value at the 2 a level. 
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Fig. 6. — Distributions (histograms) of the intermediate scale disk confidence levels. The top row shows the results for disk set A, the 
bottom for disk set B. The columns show, left to right, the two-, three-, and four-point function results. The gray bands indicate 1 and 
2 cr confidence regions, computed from 5000 simulations. The solid line indicates the results from the foreground corrected map, and the 
dotted line shows the results from the original co-added map. 



TABLE 2 

Intermediate scale results 




TABLE 3 

Small scale results 




Correlation function Both sets Disk set A 


Disk set B 


Both sets Disk set A 


Disk set B 




No foreground correction 




No foreground correction 




Two-point 0.829 0.322 
Three-point function 0.955 0.805 
Four-point function 0.941 0.917 


0.944 
0.981 
0.722 


Three-point function 0.725 0.600 
Four-point function 0.325 0.507 


0.554 
0.360 


Foreground correction by external templates 


Foreground correction by external templates 


Three-point function 0.330 0.084 


0.214 


Two-point function 0.756 0.189 


0.901 


Four-point function 0.177 0.639 


0.317 



Three-point function 0.816 0.527 0.938 
Four-point function 0.683 0.674 0.330 

Note. — Results from the intermediate scale full-sky 
tests. The numbers indicate the fraction of simulated 
realizations with value lower than that for the co- 
added WMAP map. The upper half shows the results be- 
fore correcting for foregrounds, and the lower half shows 
the results after applying foreground corrections. 



Note. — Results from the small scale full-sky tests. 
The numbers indicate the fraction of simulated realiza- 
tions with x^ value lower than that for the co-added 
WMAP map. The upper half shows the results before 
correcting for foregrounds, and the lower half shows the 
results after applying foreground corrections. 



From these considerations, it seems likely that the sim- 
ple foregromidcorrectwn method by templates discussed 
bv lBennett et al.1 l)2003bD leaves significant residuals near 



the galactic plane. Indeed, this should not be surprising 
since the input synchrotron and free-free templates do 
not contain power on the small angular scales probed by 
the Q-, and W^-band maps, and the template fitting 
method itself does not admit spectral variations on the 



sky, while it is likely that such variations are seen close 
to the galactic plane. 

We also compute the full-sky, disk-averaged correla- 
tion function for the small-scale functions, and the results 
from this analysis are shown in Table 01 Here we see 
that the model is comfortably accepted on these scales, 
and the effect of the foreground residuals discussed above 
is diluted by the additional sky coverage. 

Finally, we make one connection to a previously re- 
ported detection of non-Gaussianity ijVielva et alJl200'^ 
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Fig. 7. — The full-sky intermediate-scale correlation functions. The left-hand column shows the functions directly as measured from the 
union of the two disk sets, while the normalized functions (see equation |3 are shown in the other three columns. The gray bands indicate 
1, 2, and 3cr bands, as computed from simulations. The dotted line corresponds to the foreground corrected map and the solid line to the 
uncorrected map. 



TABLE 4 
Galactic plane results 



Scales 


Two-point 


Three-point 


Four-point 


Intermediate 


0.533 


0.826 


0.975 


Small 




0.975 


0.998 



Note. — Results from tests of the correlation 
functions computed over the disks near the galactic 
plane (disks A34-54). 



ICruz et al.ll20nl : A very cold spot was found at galactic 
coordinates {b = —57°, I = 207°) using wavelet statistics, 
this corresponds to disk number B73 (see FigureEJ in our 
partitioning of the sky. This particular disk has a three- 
point function that is high at the 2 a level on both 
intermediate and small scales, insignificant by itself, yet 
perhaps inte r esting when taken in combination with the 
iVielva et all l|2004(i detection. 

6. CONCLUSIONS 



We have computed the two-, three-, and four-point 
correlation functions from the first-year WMAP data 
sets, and find interesting effects on several angular scales. 
On the very largest scales an asymmetric distribution of 
power is observed in both the two-, three- and four-point 
functions, in that the fluctuations on the southern eclip- 
tic (and galactic) hemisphere are significantly stronger 
than on the northern ecliptic (and galactic) hemisphere. 
In order to study this effect more closely, we computed 
the correlation functions from each frequency band sep- 
arately, and for two different sky cuts, and found that 
the asymmetry is present in any of the bands, and inde- 
pendent of the particular region definition. This argues 
against a foreground based explanation for this effect. 

Next, we computed the two-point correlation functions 
from a set of difference maps, and detected excess corre- 
lations in the data among the PF-band differencing as- 
semblies, which are not accounted for in the detailed 
simulation pipeline used by the WMAP team. While 
this effect could potentially pose a serious problem for 
the upcoming polarization data, its absolute amplitude 
is very small compared to the temperature anisotropy 
amplitude, and it is therefore highly unlikely to cause 
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Fig. 8. — Results from the small-scale correlation function analysis, 
values around the galactic plane in disk set A. 



any problems for results based on the first-year WMAP 
temperature data. 

We then computed the correlation functions on small 
(< 72') and intermediate (< 5°) scales, and found that 
the agreement with the Gaussian model is generally good 
in these cases. Although a pattern consistent with the 
large-scale asymmetry discussed earlier is visible in the 
intermediate-scale three-point correlation function, it is 
difficult to assess the significance of this pattern. It 
should be regarded more as supportive evidence to the 
large-scale results, than as a conclusive result on its own. 
Overall, the Gaussian model is accepted by the interme- 
diate scale A^-point correlation functions. 

On small scales we detect residual foregrounds near 
the galactic plane roughly at the 2.5a level. However, 
such residuals may be seen by eye in the actual maps, 
and this is therefore not a surprising result. Except for 
this residual foreground detection the model is accepted 
by iV-point correlation functions on the smallest scales 
probed in this paper. 

As seen from the analyses presented in this paper, real- 
space based statistics, such as the A^-point correlation 
function, have a clear value with respect to control of 
systematics. For cosmological purposes, harmonic-space 
statistics (e.g., the angular power spectrum, and the bi- 
spectrum) are usually the preferred tools since they gen- 
erally have a simpler physical interpretation than their 
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real-space counterparts. However, systematics are often 
localized in real space rather than in harmonic space 
(e.g., foregrounds are highly localized in space; 1// noise 
leads to stripes along the scan directions; cross-talk be- 
tween detectors leads to noise correlations at some given 
scale), and real-space statistics can therefore often be 
more powerful for detecting their presence. The results 
presented in this paper are clear demonstrations of this 
fact. 
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Fig. 9. — Distributions (histograms) of the small scale disk confidence levels. The top row shows the results for disk set A, the bottom 
for disk set B. The columns show the three- and four-point function results, respectively. Gray bands indicate 1 and 2 a confidence regions, 
computed from 5000 simulations, the solid lines indicate the results from the foreground corrected map, and the dotted lines show the 
results from the raw co-added map. 
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